clc
clear
close all
%% Pulse shape & Variable initialization
MAIN             = MainFunctions;
pulse            = 4;      % 1 -> lorentzian pulse
                           % 2 -> GMSK pulse BT = 0.3
                           % 3 -> LRC pulse
                           % 4 -> LREC pulse
pulse_length     = 1;      % 1  -> Full response
                           % >1 -> Partial response
os               = 8;    % Over sampling frequency
Ts               = 1/os;   % Sampling Time
M_ary            = 2^1;    % M_ary symbols used 2 -> Binary
modulation_index = 0.5;    % Modulation index
width            = 0.7;    % This variable is used for Lorentzian Pulse only. (Not be used for pulse > 1)
dmin             = 1.78;   % GMSK BT=0.3
% 'decision_delay' Decide after how many observation symbols
% the Viterbi receiver can make a decision on the detected bit.
decision_delay   = 50;
%-------------- Modulated data ----------------------
snr           = 30:30;
BER_all       = zeros(1,length(snr));
%% frequency pulse
[g_t,q_t]     = MAIN.CREATECPMPULSE(pulse,pulse_length,width,os,0); % Function return the CPM pulse and phase.
                                                         % g_t = g(t) is the CPM pulse shape.
                                                         % q_t -> is the phase, integral of g_t.                                                      
m  = 1:M_ary;
if(pulse~=1)
    A_m    = 2*m-1-M_ary;       % Different modulation level
else
    modulation_index      = 2*modulation_index;
    A_m    =  0:M_ary-1;        % Different modulation level
end
%-------------------------------------------------------------------------------------------------------------

%% ---------------------------------------------- Receiver------------------------------------------------------
%---------------------------- states -- Branches -----------------------------------
[States_Number,Branche_Number,phase_states] = MAIN.states_Branches(modulation_index,pulse_length,M_ary);
% [States_Number,Branche_Number,phase_states] = MAIN.states_Branches_rimoldi(modulation_index,pulse_length,M_ary);%加倾斜相位之后的状态分支

[Branches,states_sort] = MAIN.states_transition(States_Number,Branche_Number,phase_states,pulse_length,M_ary,A_m,modulation_index);
% [Branches,states_sort] = MAIN.states_transition_rimoldi(States_Number,Branche_Number,phase_states,pulse_length,M_ary,A_m,modulation_index);%加倾斜相位之后的状态转移
%------------------------------- Branche_metric ----------------------------------
if(pulse_length>1)
    Branche_metrics = zeros(Branche_Number,length(((pulse_length-1)*os:pulse_length*os)));
else
    Branche_metrics = zeros(Branche_Number,length(1:pulse_length*os+1));
end
for i = 1:Branche_Number
    alpha                = Branches(i,2:end);
    Branche_metrics(i,:) = MAIN.compute_Branche_metric(os, modulation_index, pulse_length, alpha, g_t);%计算8个分支对应的分支矩阵
end

    
%-------------theoretical BER calculation-----------------------
EbNodB_vect_lin = 10.^(snr/10);
x               = sqrt(dmin.*EbNodB_vect_lin);
Q               = 0.5.*erfc(x/sqrt(2));
Pe              = Q;
%---------------------------------------------------------------
frame_max_length = 300./Pe*log2(M_ary);
BER_all    = zeros(1,length(snr));
tic
NBits_limit = [10^4,10^5]; % min, max of number of bits in the whole frame
% Itterate over all the snr
for snr_no = 1:length(snr)
    NBits      = frame_max_length(snr_no); % Number of bits (number of bits the whole frame)
    if(NBits>NBits_limit(2))
        NBits = NBits_limit(2);
    end
    if(NBits<NBits_limit(1))
        NBits = NBits_limit(1);
    end
    pack       = 10^4;                    
    Nbits_div  = round(NBits/pack); % Number of packs in a frame (each = to 10^4)
    % For each snr iterations, iterate over all Nbits_div, where
    % Nbits_div is the number of packs.
    for I = 1:Nbits_div
        Nbits       = log2(M_ary)*floor(round(NBits/Nbits_div)/log2(M_ary));
        rand('seed', 16);
        bits  = randi([0 1],Nbits,1).';
        bits_m = MAIN.Mapp(bits, M_ary, pulse);
        Nbits  = length(bits_m);               % Number of bits change after applying the mapping (only for M-ary>2).
        %% --------------------------- Modualtion -----------------------------------
        %----------------- M odulated signal (information dependent)------------------
        [received_signal, tx_signal] = MAIN.Modualtion(bits_m,os,Ts,Nbits,g_t,modulation_index,pulse_length,snr_no,M_ary,snr,q_t);

        %------------------ Viterbi Decoding ------------------------------ 
        [DetectedBits] = MAIN.viterbiDecoding(Nbits,decision_delay,States_Number,pulse_length,Branche_Number,Branches,received_signal,os,Branche_metrics,Ts,states_sort,M_ary);

        if(M_ary==2)
            BER(snr_no)       = nnz(bits_m(2+pulse_length:end-decision_delay-2-pulse_length)-DetectedBits(1+pulse_length:end-3-pulse_length));
        elseif(M_ary==4)
            BER(snr_no)       = nnz(bits_m(2+3*pulse_length:end-decision_delay-3-3*pulse_length)-DetectedBits(1+3*pulse_length:end-4-3*pulse_length));
        elseif(M_ary==8)
            BER(snr_no)       = nnz(bits_m(2+3*pulse_length:end-decision_delay-3-3*pulse_length)-DetectedBits(1+3*pulse_length:end-4-3*pulse_length));
        end
        
        BER_all(snr_no) = BER_all(snr_no)  + BER(snr_no);
        sprintf('BER_all for division number : %d / %d for the snr %d',I,Nbits_div,snr(snr_no))
        if(I==Nbits_div) % at the last frame
            if(M_ary==2)
                BER_all(snr_no) = BER_all(snr_no)/(length(bits_m(2+pulse_length:end-decision_delay-2-pulse_length))*Nbits_div);
            elseif (M_ary==4)
                BER_all(snr_no) = BER_all(snr_no)/(length(bits_m(2+3*pulse_length:end-decision_delay-3-3*pulse_length))*Nbits_div);
            elseif (M_ary==8)
                BER_all(snr_no) = BER_all(snr_no)/(length(bits_m(2+3*pulse_length:end-decision_delay-3-3*pulse_length))*Nbits_div);
            end
            sprintf('BER is %f',BER_all(snr_no))
        end
        
    end
end
toc
%----------------------------------Plots----------------------------------------------
figure(1)
semilogy(snr, BER_all, 'linewidth', 3);hold on;
semilogy(snr,Pe,'linewidth', 3);
xlabel('E_b/N_0 (dB)'); 
ylabel('BER');
legend('Simulation', 'Bound');
grid on;
